$NOFLOATCALLS
$NODEBUG
$STORAGE:2
c*****************************************************************************
c	NOTE: Some of the variables are not the same as in the article
c	      For example hadj = v, est=e(i+1), hfact=s, hlim = vlim.
c*****************************************************************************
c
c
	subroutine runkut(xa,ya,xb,neqn,tola,tolr,hstart,w,
     &			  imeth,ierror,icom,func)
c
c****	This is simply a front-end subroutine to split up the work
c****	array specified by the user into smalller arrays
c
	implicit double precision (a-h,k,o-z)
	external func
	dimension icom(4),w(*)
	common /epsil/eps

	call solver(xa,ya,xb,neqn,tola,tolr,hstart,ierror,imeth,
     &		    icom(1),icom(2),icom(3),icom(4),w(1),w(1+neqn),
     &		    w(1+neqn*2),w(1+neqn*3),w(1+neqn*4),func)
	return
	end
c*****************************************************************************
c
c
	subroutine solver(xa,ya,xb,neqn,tola,tolr,hstart,ierror,imeth,
     &			  ist,ichk,idef,iround,kt,est,yold,ynew,w,func)

	implicit double precision (a-h,k,o-z)
	logical hflag
	external func
	dimension ya(neqn),w(13,neqn),kt(neqn)
	dimension est(neqn),yold(neqn),ynew(neqn)
	common /cons/hmin,hmaxl,hfact,hlim,power
	common /epsil/eps

c****	If TOLA is set to zero, then a relative error test
c****	If TOLR is set to zero, then an absolute error test
c****	If neither are set to zero, then a mixed error test
c****	If this is the first call to RUNKUT calculate the machine
c****	epsilon and work out the constants depending on method used.
c****	To ensure that the results will always be "safe", the value
c****	EPS used by RUNKUT is actually 20 times the machine epsilon.

	if(ist.EQ.0) then
	  call caleps
	  call const(ist,imeth,idef)
	  eps=eps*20.d0
	end if

c****	Check if the sum of TOLA and TOLR is greater than 10*EPS
c****	If not, set TOLR to one million times EPS.

	ierror=0
	if(ichk.GT.0)then
	  if(tola+tolr.LT.10.d0*eps)then
	    ierror=1
	    tolr=1.d06*eps
	  end if
	end if
	iround=0

	isig=dint((xb-xa)/dabs(xb-xa))
	hold=hstart
	x=xa
	do 10 i=1,neqn
	   yold(i)=ya(i)
10	continue

c****	Set HMAX to the maximum value HMAXL set in CONST unless this
c****	is greater than the interval width, in which case HMAX is
c****	equal to the interval width.

	hmax=dmin1(hmaxl,dabs(xa-xb))

c****	call the runge-kutta solver using the current step size
c****	to predict y(n+1)

15	call rkp(x,ya,hold,neqn,imeth,kt,est,yold,ynew,w,func)

c****	calculate step size adjustment factor HADJ
c****	then calculate HNEW using the smallest value of HADJ

	hadj=9999.d0
	do 20 i=1,neqn
	   absest=dabs(est(i))
	   if(absest.EQ.0)then
	     hadjl=hlim
	   else
	     if(yold(i).EQ.0)then
	       hadjl=((tola+tolr*tolr)/absest)**power
	     else
	       hadjl=((tola+tolr*dabs(yold(i)))/absest)**power
	     end if
	   end if
	   if(hadjl.LT.hadj) hadj=hadjl
20	continue

c****	adjust the step size to HNEW using the calculated value of HADJ
c****	unless HADJ is greater than HLIM, in which case use HLIM
c****	If HNEW is larger than HMAX, choose HMAX as the new step size
c****	If HNEW is less than HOLD/10. keep HNEW as HOLD/10 to avoid
c****	excessively large swings in the step size

	hold1=dabs(hold)
	hnew = dmax1(hold1/10.,
     &		     dmin1(hfact*hold1*(dmin1(hadj,hlim)),hmax))

c****	Check if HOLD is large enough compared to EPS to avoid
c****	severe round-off errors. If HNEW is less than HMIN, the
c****	problem has got too stiff or is discontinous-exit saving
c****	the last points calculated.  If the last step was sucessful
c****	let YOLD=YNEW and calculate YNEW using new step size. If it
c****	was unsuccessful, recalculate YNEW using reduced step size.
c****	If the HOLD was a reduced step size, then restrict HNEW to HOLD
c****	If XB will be reached or exceeded,  exit after calculating
c****	Y at XB.

	if(dabs(x).GT.0) then
	  if((hnew/(dabs(x)*18.d0)).LE.eps) iround=1
	end if
	if(hnew.GE.hmin) then
	  hnew=isig*hnew
	  if(hadj.GE.1.d0) then
	    if(isig*(x+hold-xb).LT.0.d0) then
	      x=x+hold
	      if(.not.hflag) hnew=hold
	      hflag=.true.
	      do 30 i=1,neqn
		 yold(i)=ynew(i)
30	      continue
	      hold=hnew
	      goto 15
	    else
	      hstart=hnew
	      hold=xb-x
	      call rkp(x,ya,hold,neqn,imeth,kt,est,yold,ynew,w,func)
	      do 50 i=1,neqn
		 ya(i)=ynew(i)
50	      continue
	    end if
	  else
	    hflag=.false.
	    hold=hnew
	    goto 15
	  end if
	else
	  hnew=isig*hnew
	  ierror=ierror+2
	  xb=x
	  do 60 i=1,neqn
	     ya(i)=yold(i)
 60	  continue
	end if
	return
	end
